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SUMMARY 

Birkhoff's  normalizing  canonical  transformation  at 
an  equilibrium  of  elliptic  type  with  no  internal  resonance 
can  be  built  explicitly  and  recursively,  without  partial 
inversions  or  substitutions,  by  means  of  Lie  transforms. 

Invariant  sections  and  ordinary  families  of  periodic 
orbits  for  truncated  normalized  systems  are  analyzed  in 
detail. 


Introduction 


We  consider  a  conservative  dynamical  system  with  m  degrees  of 
freedom  whose  equations  of  motion  we  write  in  vector  form 


!?*• 


(1) 


in  which  q  and  p  indicate  vectors  with  components  q^  and  p^ 

(1  j  ^  m)  respectively.  We  assume  that  the  origin  (p  ■  0 ,  q  *  0) 
in  the  phase  space  is  an  equilibrium  point,  and  that  the  Hamiltonian 
3f(q,p)  is  an  analytic  function  of  the  state  variables  in  a  neighbor¬ 
hood  of  the  equilibrium.  Hence  Taylor's  expansion  of  3i  in  that 
neighborhood  is  of  the  form 


3f(p,q)  “  ^  ~T  5f_(p,q)  , 

n>0  nI  n 


(2) 


where,  for  each  n  0,  is  a  homogeneous  polynomial  of  degree 

n  +  2  in  the  state  variables.  Indeed  the  constant  term  in  the  develop¬ 
ment  of  31  around  the  equilibrium  does  not  contribute  to  the  equations 
(1),  and  since  (q  -  0,  p  =  0)  is  a  solution  of  (1),  the  first  degree 
terms  in  (2)  must  vanish.  Thus  the  power  series  for  31  begins  with 
quadratic  terms. 

We  suppose  that  the  series  (2)  has  only  real  coefficients,  in  particular 
that  the  quadratic  form  31  q  is  real.  Let  be  the  matrix  of  the  linear 
Hamiltonian  system  derived  from  3{q.  If  X  is  an  eigenvalue  of  ,  then 
its  complex  conjugate  X  and  its  inverse  are  also  eigenvalues  of 


Of.  Tims  assuming  that  all  eigenvalues  of  0\  are  distinct  and 
purely  imaginary,  we  can  represent  them  by  the  sequence 


iVj.iX.,,...,!!^,— i\^,— i,\ 


■l\  , 

in 


(3) 


where  \  are  non  zero  real  numbers.  Under  these  conditions, 

J  J.  m 

there  exists  a  canonical  transformation 


qj "  JLbA  “s  *k  *  8jA  sln  ♦k|> 

pj  ■  cos  *k +  sjA 81,1  \| 


(4) 


from  the  original  state  variables  (q,p)  to  angle  coordinates  $  =  (^) 

and  action  momenta  1  =  (1.  )  such  that  the  quadratic  form  converts 

K  0 

into 


•  .  r  V 


(5) 


l^kMn 


Under  the  mapping  (4),  the  homogeneous  component  transforms  into 

a  homogeneous  polynomial  of  degree  n  +  2  in  the  state  variables 
its  coefficients  being  trigonometric  sums  in  the  arguments  | ^  witli 
real  coefficients.  For  each  n  2.  the  component  3f  (1,^)  lias  the 
d'Alenbcvt  ahaiviatcinatic.  It  means  that,  for  a  term  like 

(y/^l)  l(^)  ‘  (An)  cos(kl*l  +  k2*2  +  +  kn,0 


to  enter  the  homogeneous  component  of  degree  p  +  p +  •  •  •  +  p  ,  it 
is  necessary  that,  for  1  *-  J  v  m. 
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I k j |  ^  p.  and  k^  p(  (mod.  2);  (6) 

the  same  condition  is  lo  be  satisfied  by  any  sine  term. 

Birkhol'f's  normalization  consists  in  finding  for  a  given  N  ^  1 
a  canonical  transformation  0:(l,-|>)  — >  (I*,|*)  such  that  the  converted 
Hamiltonian  decomposes  into  the  sum 

ard*.**)  -ovci*,-)  +  (7) 


where  cX  is  a  polynomial  with  real  coefficients  in  the  transformed 
actions  1*  =  (I*)  of  degree  at  most  equal  to  (N  +  2)/2,  and 

K 

is  a  series  in  the  action  amplitudes  vl*.  beginning  wi tli  terms  of  at 
least  degree  (N  +  3).  the  coefficients  being  periodic  functions  of 


the  angles 


"k- 


Birkhoff  t.1927)  has  shown  that  the  normalized  part  rV(I*,-)  of 
the  Hamiltonian  is  invariant  with  respect  to  the  various  canonical 
t rans format  ions  leading  to  it.  Hence  it  yields  dynamical  characteristics 
of  the  flow  of  trajectories  close  to  the  equilibrium. 


By  halting  the  sequence  of  normalizing  transformations  after  a 
finite  number  N  of  steps,  we  obtain  an  asymptotic  representation  of 
the  phase  portrait  at  the  equilibrium  as  the  actions  I*  — >  0.  Practical¬ 
ly  the  original  system  as  described  by  .‘Tt  is  replaced  by  its  asymptotic 
approximation  truncated  alter  order  N  +  2  as  represented  by  c'V.  But, 
as  the  angles  are  ignorable  in  rV,  the  dynamical  system  that  it 

describes  is  integral!  le,  and  its  structure  in  phase  space  is  the  product 
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of  an  m-dimensional  torus  by  the  m-dimensional  phase  space  (Arnold 

and  Avez  1967).  Now  most  of  the  invariant  tori  of  quasi-periodic 

motions  in  the  flow  of  trajectories  determined  by  oV  will  survive 

(N) 

the  perturbations  induced  by  the  remainder  ' .  This  may  explain 
why  (although  Birkhoff’s  normalization  is  only  exceptionally  convergent 
as  N  —>“>),  it  proves  useful  in  revealing  the  dynamical  structures 
around  the  equilibrium  for  the  full  system 


1.  Normalization  by  Lie  transforms 

It  is  assumed  that  the  basic  frequencies  X..  ,X X  are 

1  L  m 

rationally  independent  of  each  other. 

Traditionally  the  normalizing  transformation  is  defined  from  a 
generating  function 

W(I*,0)  -  I  I**.  +  I  W  (I*,*).  (8) 

1 <k<m  K  *  n> 1  n 

by  the  implicit  equations 


*k  ’  *k  +  31*  Wn’ 

n>l  k 


I  -  I*  +  V  —  W  . 

k  k  n>l  ^k  n 


(1  <  k  <  m) 


(1  <  k  <  ra) 


(9) 


Conversion  is  then  made  to  the  complex  variables 


Ck  "  y/\  6 


"k  ■  -1  A 


(1  <  k  <  m)  (10) 


so  that  the  operations  of  normalization  are  performed  entirely  within 
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the  algebra  of  polynomials  in  2ra  letters  over  the  complex  numbers. 

They  have  been  programmed  for  automatic  processing  by  computers  for 
systems  with  two  degrees  of  freedom:  programs  originally  set  up  to 
deal  with  resonance  cases  (Gustavson  1966)  have  been  modified  to 
operate  in  the  regular  case  (Chai  and  Kass  1966). 

The  use  of  a  generating  function  in  order  to  build  a  normalizing 
canonical  transformation  is  reminiscent  of  a  well  known  procedure  in 
the  theory  of  General  Perturbations  (Poincard  1893)  that  is  usually 
referred  to  as  Von  Zeipel's  method.  But  its  shortcomings  have  recently 
prompted  the  formulation  of  a  more  direct  algorithm  based  On  Lie  trans¬ 
forms  (Deprit  1969).  This  new  formalism  simplifies  to  a  large  extent 
Birkhoff’s  normalization.  The  present  application  will  in  turn  check 
the  general  procedure.  For  we  had  previously  normalized  the  Hamiltonian 
of  the  restricted  problem  of  three  bodies  at  (Deprit  et  at.  1967), 

and  we  shall  find  that  the  algorithm  about  to  be  described  indeed  restores 
the  already  established  normalized  Hamiltonian. 

The  basic  step  is  to  construct  a  sequence  of  functions 
WjU*,;*)  (l^j^N) 

depending  not  on  mixed  variable®  but  on  both  new  angles  and  new  actions. 
They  will  be  obtained  in  the  course  of  transforming  the  Hamiltonian.  For 
reasons  of  clarity,  let  the  exposition  refer  to  the  diagram  in  Fig.  1 
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Fig.  1.  -  Recursive  normalization  of  the  Hamiltonian. 


The  functions  entering  the  diagonal  5£q ,  3f^,3l 2 » • • .  are  the  homo¬ 
geneous  components  in  the  initial  series  (2)  as  expressed  by  means  of 
the  primitive  angles  and  actions  ^  and  Ij  (1  _<  j  <_  m) ;  the  functions 
entering  the  diagonal  are  their  normalized  forms. 

The  elements  in  Table  I  are  to  be  built  recursively  according  to  the 

law 


«K(p)  -*X(p_1)  +  V  /‘WP'1)*W  \ 

Vl  n<7<a  J+V 


where 


0<J_<q 

j  is  the  binomial  coefficient 


(3) 


q(q-l)  •••  (q-i-H) 


1.2. 


j 


(ID 


and  l  is  the  Poisson  bracket 

'  q-j  j+1' 

/n/( P"D.  y  \  *  \  /  ■ «/(P”l)  — —  U  -  — 9£^P_ 

Pq-j  ’  Vl/  K^m^k  q"j  3Ik  j+1  3Ik  ^ 


-1)  3 


W. 


j+lj 


Extending  a  proposition  established  by  Brown  and  Shook  (1933) ,  we  can 
prove  that,  if  f  is  a  homogeneous  d'Alembert  series  of  degree  p  with 


respect  to  the  variables  v/l*,  yi£ , . . . ,  . . .  ,<(>*  and  g 

is  a  homogeneous  d'Alembert  series  of  degree  q  with  respect  to  the 

same  variables,  then  the  Poisson  bracket  (f;g)  is  also  a  d'Alembert 

series  with  respect  to  these  variables,  but  of  degree  p  +  q  -  2.  The 

fact  is  that,  while  the  separate  partial  derivatives  with  respect  to 

the  actions  I* I*,..., I*  are  not  d'Alembert  series,  in  each  Jacobian 
12  m 

a(f.g)  .  1  3(f.«) 

2  /if 

the  divisor  yij*  disappears,  and  by  exact  cancellation  of  all  terms 
not  having  the  d'Alembert  characteristic,  the  Jacobian  of  f  and  g 
with  respect  to  <f>*  and  yT*~  turns  out  to  be  a  homogeneous  polynomial 
in  the  variables  y  I*,  yi*\. . . ,  yT*"  of  degree  p  +  q  -  2  with  the 
d'Alembert  characteristic. 


This  proposition  assures  the  consistency  of  the  scheme  we  propose 
here;  it  is  useful  as  a  check  on  the  various  parts  of  a  program  to  be 
implemented  by  computer. 


As  we  shall  see,  the  differential  operator 


DO) 


V 


1  7  k  3$,* 

l<k<m  k 


will  play  an  essential  role.  It  is  a  linear  mapping  into  itself  of  the 

algebra  <3f  <5f(I*,£*)  of  d'Alembert  series  with  respect  to  yf*,  yi*, . . . , 

yT*,  . 4>*.  More  precisely  it  is  a  linear  mapping  into  itself 

v  m  12  m 

of  the  vector  space  (I*  ,4>*)  of  homogeneous  polynomials  of 
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degree  p  in  yi* ,  yi*, . . . , yi*  having  the  d'Alembert  characteristic 
and,  under  the  condition  that  the  frequencies  X  ■  (X^)  are  rationally 
independent  of  each  other,  the  kernel  of  this  restriction  is  the  vector 
space  h  (1*)  of  homogeneous  polynomials  in  yi*,  yi*, . . . ,  yi* 

of  degree  p  with  real  coefficients.  Notice  also  that 


D(X)cos/  £ 

'lf_k<m 

D(X)sin/ 

l<k<m 


(12) 


which  makes  the  inverse  operator 


(D(X))"1 


a  trivial  operation  over  &f. 


Now  let  us  see  how  the  normalization  is  carried  out  recursively. 

In  view  of  (11) , 

^oX)  wL),  (13) 

with  the  understanding  that  in  J(q  and  the  letters  I  and  <J> 

have  been  replaced  by  I*  and  <f>*  respectively.  We  collect  in 

the  terms  of  3^  that  do  not  depend  explicitly  on  the  angle  coordinates 

(3) 

4)*.  Of  course,  since  3^  is  an  element  of  Or  ' ,  there  are  no  such 
terms,  and  we  find  that  3^^  must  be  put  equal  to  zero.  On  the  other 
hand 


Wx)  -  -  D(X)Wr 

so  that  the  definition  (13)  transforms  into  the  linear  partial  differential 
equation 

D(X)WX  -3fr  (14) 
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A  particular  solution 


W, 


d(IFwi 


of  (14)  is  constructed  according  to  the  rules  (12),  and  thus  it  is  an 
(2) 

element  of  ' .  We  neglect  to  complete  it  by  arbitrary  elements  of 
the  kernel  of  D(A);  their  addition  would  indeed  deny  the  d’Alembert 
characteristic  to  W^. 


Suppose  that  the  first  N  rows  (N  >  1)  of  the  triangle  in  Fig.  1 
1  have  been  computed,  that  is  to  say  the  Hamiltonian  has  bean  normalized 
up  to  order  N  +  1  and  the  first  N  -  1  generators  W^ ,W2» . . .  ,WN_^  of 
the  canonical  transformation  have  been  computed.  At  this  step,  it  is 
convenient  to  put  WM  «  0  for  the  time  being,  and  to  compute  the  elements 
of  the  (N  +  l)-th  row  of  the  triangle  in  Fig.  1  by  the  formula  (11). 

***  (n) 

Denoting  them  by  3i  ,  we  obviously  have  that,  for  p  >_  1  and  p  +  q  -  N, 


*qP)  ”5qP>  +  «0!  V 


<£?(p) 


Notice  that  each  element  belongs  to  the  vector  space  Cf 


(15) 

(N+2) 


Making  p  =  N  and  q  =  0  in  (15),  we  find  the  relation 


v  -ar™.  (16) 

(N)  ^(N) 

Once  more  we  collect  in  the  terms  of  3^  which  do  not  depend 

(N) 

explicitly  on  the  angle  coordinates  <J>*.  If  N  is  odd,  then 

(N) 

turns  out  to  be  zero;  but  if  N  =  2M  is  even,  is  a  homogeneous 

polynomial  of  degree  M  in  the  letters  I*, I*,..., I*  with  real  coef¬ 


ficients.  Let  us  put 


zzrrr. 
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j,(n)  ,as> 


(17) 


so  that  we  can  write  the  identity  (16)  as  the  linear  partial  differential 
equation 

D(X)Wn  =  JP(N\  (18) 


of  which  a  solution  is  found  in  accordance  with  the  rules  (12)  : 


=  -JL  *00 

N  D(A)  *  ‘ 


(19) 


Observe  that  the  particular  solution  (19)  belongs  to  the  vector  space 
<^(N+2) .  if  n  is  odd,  completing  with  elements  of  the  kernel  of 

D(\)  would  result  in  depriving  the  completed  solution  of  its  d'Alembert 
characteristic.  But,  if  N  =  2M  is  even,  that  characteristic  would  not 
disappear  by  adding  to  WN  a  homogeneous  polynomial  of  degree  M  in 
the  letters  I*, I*,..., I*  with  real  coefficients.  Yet,  as  follows  from 
Birkhoff's  proof  of  invariance  with  respect  to  the  group  of  normalizing 
canonical  mappings,  this  complementary  term  is  not  going  to  contribute 
subsequently  to  the  normalized  Hamiltonian.  Therefore  we  decide  to  omit 
it  systematically. 

Having  determined  W^,  we  complete  the  elements  (15)  as  follows 

^qP)  =^qP)  -'*(N)  (p  >  1,  q  >  0,  p  +  q  =  N)  (20) 

so  as  to  be  ready  to  extend  the  normalization,  if  necessary,  beyond  the 
order  N  +  2. 


So  far  we  have  determined  the  components  of  the  normalized 
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Hamiltonian,  and  the  generators  . .  ,WN)  of  the  Lie  transform. 

In  order  to  construct  formally  the  normalizing  canonical  transformation, 
we  have  to  compute  the  Lie  transforms  of  the  original  angle  and  action 
variables  (<f>,I).  To  this  effect,  we  apply  the  usual  triangular  algorithm 
(Deprit  1969)  which  converts  series  like 

*>  =<£,  tp  E  0  for  n  >  1 

into  functions  of  the  normalizing  angle  and  action  variables  ( <p* ,  I*) . 

We  proceed  in  the  same  manner  to  express  the  original  state  variables 
(q,p)  as  given  by  (4),  if  we  need  them  expressed  in  terms  of  (<p* , I*) . 

2.  Normalized  systems  with  two  degrees  of  freedom 

From  now  on  we  restrict  o  -rselves  to  dynamical  systems  with  two 
degrees  of  freedom 

As  we  just  described  it,  Birkhoff’s  normalization  has  now  become 
a  relatively  easy  operation  to  be  implemented  automatically  by  computer, 
and  the  analysis  of  an  elliptic  equilibrium  can  avail  itself  freely  of 
this  routine  technique.  Hence,  from  the  methodological  standpoint,  the 
emphasis  shifts  on  the  problem  of  how  to  obtain  local  information  about 
the  phase  space  at  the  equilibrium  from  the  normalized  Hamiltonian. 

Quite  naturally  we  are  led  to  contrast  two  kinds  of  dynamical 
systems.  On  one  hand  there  is  the  original  system  described  by  the 
full  Hamiltonian  which,  after  transformation,  takes  the  form 
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(7);  on  the  other  hand,  at  a  given  order  N  +  2,  there  is  the  integrable 
system  described  by  the  principal  part  oV(I*,-).  Before  we  indicate  how 
the  integrable  system  is  a  deformed,  yet  informative,  version  of  the  full 
system,  let  us  see  how  we  can  study  in  depth  the  principal  part. 

The  canonical  transformation  0* :  (q1  ,q2  ,p^ ,p2)  — >  (<£*,<(>*, I*, I*) 
from  the  original  state  variables  to  the  normalizing  angle  and  action 
variables  is  singular  along  the  phase  manifolds  I*  5  I*^CJi»q2 ,P1,P2^  *  ® 
and  I*  =  I^(q1,q2,P1,P2)  -  0. 

But,  for  the  normalized  system,  i.e.  the  dynamical  system  represented 
by  the  integrable  Hamiltonian  (H ,  the  equilibrium  configuration  is  defined 
by  I*  =  I*  =  0.  Moreover,  the  frequencies  X^  and  *2  being  rationally 
independent,  there  emanate  from  the  equilibrium  two  natural  families  of 
periodic  orbits.  The  f airily  <7-^  associated  with  the  frequency  X^  lies 
on  the  manifold  I*  =  0,  and  conversely  any  point  of  the  integral  mani¬ 
fold  I*  =  0  belongs  to  one,  and  only  one,  periodic  orbit  in  the  family 
£7^.  Likewise  the  family  $F2  that  Liapunov's  theorem  associates  with 
the  frequency  X2  lies  on  the  manifold  I*  *  0,  and  it  is  coextensive 
to  that  set. 

It  is  easy  to  show  (Deprit  et  at ,  1967)  that  the  nontrivial  characteris¬ 
tic  exponents  for  the  periodic  orbits  belonging  to  £7^  and  are 

respectively  the  series 


t 
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±i(X^  +  a^I*  +  a2I*2  +  •••) 
±i(X2  +  b  I*  +  b2I*2  +•••) 


in  the  powers  of  I*  or  I*.  The  independence  of  X^  and  X2  over 
the  rationals  implies  of  course  that  neither  X^  nor  X2  is  equal  to 
zero.  Consequently  the  nontrivial  multipliers  of  the  periodic  orbits 
either  in  0^  or  in  0^  are  not  equal  to  one;  in  other  words  for  these 
orbits,  the  algebraic  multiplicity  of  the  multiplier  +1  is  equal  to  2, 
which  implies  that  the  gradients  of  the  integrals  I*  and  I*  in  the 
original  phase  space  (q^ ,q2 .p^ ,P2)  at  any  point  of  the  integral  mani¬ 
fold  I*  =  0  and  I*  *  0  are  collinear. 

Outside  the  manifolds  I*  =  0  or  I*  =  0,  the  normalizing  canonical 
transformation  0*  is  regular,  i.e.  its  Jacobian  matrix  is  invertible. 

As  a  consequence,  since  the  gradients  of  I*  and  I*  in  the  phase  space 
(qi ,q2 ,Pi,P2)  form  two  rows  of  the  Jacobian  of  (0*)  \  they  must  be 
linearly  independent.  In  other  words,  along  any  solution  of  the  normalized 
system,  which  is  neither  the  equilibrium  nor  an  orbit  of  ^  or  0^ ,  the 
gradients  of  the  "adelphic"  integrals  are  not  collinear. 

A  torus  (I*  4  0,  I£  4  0)  is  filled  by  quasi-periodic  motions  having 


the  frequencies 
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O’-blts  In  that  torus  are  transformed  into  each  other  by  properly  adjusting 
the  additive  constants  in  the  phase  angles  .f*  and  4>i£.  In  composition 
with  the  normalizing  canonical  transformation  t'*,  this  translation  of 
origin  in  the  torus  constitutes  what  Whittaker  (1916)  called  the  aJt'lphis 
t  iww j'owkxt  :\»i. 

Exceptionally  the  constant  values  given  to  t he  integrals  I*  and 
1*  can  be  sucli  that  the  frequencies  are  commensurable:  there  could 
exist  a  pair  (p.  , p.})  of  relatively  prime  integers  sucli  that 

P2V1(11’12)  “  PlV2(1l*12)'  (21) 

In  this  case,  the  quasi-periodic  motions  degenerate  into  the 
families  of  periodic  orbits  mentioned  by  Whittaker  (1916).  For  such 
orbits,  the  fact  that  the  adelphic  integrals  have  linearly  independent 
gradients  implies  that  all  multipliers  are  equal  to+1;  in  each  torus, 
all  the  orbits  have  the  same  period. 

The  totality  of  motions  for  the  normalized  system  is  aptly  sum¬ 
marized  in  the  plane  of  actions  (I*, I*)  (see  Fig.  2).  The  origin 
represents  the  equilibrium;  to  each  point  along  the  I*-axis  corresponds 
one  and  only  one  singular  periodic  orbit  in  the  family  £7j ;  likewise 
to  each  point  of  the  L*-axls  corresponds  one  and  only  one  singular 
periodic  orbit  in  the  family  £7.,.  Elsewhere  in  the  plane  of  actions, 
the  correspondance  between  pairs  (1 *,!.*)  and  orbits  of  the  system  is 
any  point  (1**1*)  not  on  the  action  axes  represents 
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Fig.  2.  -  Tln<  plant'  ot  notions  lor  tin-  normalisation  (trimoatod  at 
onloi'  l-'i)  ot  tho  rest  rioted  problom  around  T,  in  tin'  system  sun- 
Jupiter. 
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the  class  of  quasi-periodic  motions  transformed  into  one  another  by 
idelphic  transformation. 

There  is  a  natural  network  of  curves  coordinatizing  the  plane 
of  actions  (I*,I£).  On  one  hand  the  energy  manifolds  appear  in  the 
diagram  as  the  one-parameter  family  of  curves  defined  by  the  equation 

eV(I*,I*)  ■  constant. 

On  the  other  hand  the  ratio  of  frequencies 

p  *  v(I*,I*)/v2(I*,I*), 

which  we  found  suggestive  to  call  the  rotation  number,  defines  another 
family  of  curves.  When  p  is  rational,  each  point  of  such  a  curve 
represents  the  family  of  ordinary  periodic  orbits  filling  the  torus 
(I*, I*).  We  shall  modify  Whittaker's  definition,  and  call  ordinary 
family  of  periodic  orbits  the  two-parameter  set  0(p2/p-^)  of  ordinary 
periodic-orbits  corresponding  to  the  curve  (20)  in  the  plane  of  actions. 
Let  us  recall  that,  along  0( Peri°d  is  a  function  of  the 
energy  alone  (Deprit  and  Henrard  1968,  pp.  62-63). 

Much  expertise  is  gained  from  studying  invariant  sections  through 
the  phase  space  of  a  normalized  system. 

The  figures  that  follow  pertain  to  the  normalization  (truncated  at 
order  14)  for  the  restricted  problem  of  three  bodies  around  the  triangular 
center  of  libration  L^  in  the  system  sun-Jupiter.  The  state  variables 
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are  the  Cartesian  coordinates  and  velocities  in  the  synodical  frame 
of  reference  located  at  L^.  The  sections  are  cut  by  the  hyperplane 
x  =  0  (the  x-axis  is  in  the  direction  of  the  line  of  syzygies 
oriented  from  sun  to  Jupiter). 

A  plane  section  of  a  three-dimensional  energy  manifold  is  a  two- 
dimensional  surface,  which  we  project  on  the  coordinate  plane  (y,y) 
by  lines  parallel  to  the  x-axis. 

One  way  of  making  this  surface  apparent  in  its  projection  is  by 
plotting  the  projection  of  some  of  its  curves.  In  the  case  of  an 
integrable  system,  the  most  informative  curves  are  obviously  the  traces 
of  the  invariant  tori  of  quasi  periodic  orbits  in  the  energy  manifold, 
in  particular  the  tori  corresponding  to  ordinary  families  of  periodic 
orbits.  Fig.  3  displays  for  a  normalized  system  sections  of  several 
ordinary  families  lying  on  the  energy  manifold  at  the  equilibrium. 

The  restricted  problem  at  for  the  system  sun-Jupiter  was  normalized 

up  to  degree  14;  the  residual  in  the  Hamiltonian  was  neglected.  Then 

the  curves  were  computed  from  the  series  for  the  state  variables 

•  • 

x,y,x,y  in  terms  of  1* ,|*, I* ,1* .  The  integral  of  energy  defines  a 
region  of  the  plane  (y,y)  that  will  not  be  covered  by  the  projection 
of  the  section  (the  area  shaded  in  Fig.  3).  On  each  curve  that  projects 
the  intersection  of  a  natural  family  of  periodic  orbits  with  the  x-plane 
passing  through  the  equilibrium,  the  points  at  which  x  0  are  joined 
by  solid  lines — they  are  located  above  the  projection  plane  (y,y);  the 


Fig.  3.  -  Projection  on  the  plane  (y,y)  of  the  section  in  the- 
energy  manifold  C  =  3  by  the  x-plane  going  through  L^. 
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points  at  which  x  0  are  joined  by  dashed  lines.  As  can  be  seen 
from  Fig.  3,  invariant  curves  may  have  intersecting  projections.  But 
the  points  having  the  same  projection  correspond  to  distinct  phase 
states:  on  one  curve,  x  is  >  0  whereas,  on  the  other,  if.  is  <0. 

In  numerical  explorations  of  invariant  curves,  it  has  been  the  rule 
for  some  time  to  plot  only  the  intersection  points  at  which  x  is  >0; 
Fig.  3  makes  it  obvious  that,  on  adopting  this  convention,  one  discards 
so  much  information  that  a  correct  interpretation  of  the  sections  could 
be  in  peril.  We  see  for  instance  that  the  locus  of  intersection  points 
with  x  >  0  separates  into  two  disconnected  arcs  located  on  both  sides 
of  L^;  the  same  is  true  for  the  curve  of  intersection  points  with  x  <  0. 
However,  in  each  quadrant,  the  arc  corresponding  to  x  >  0  falls  in 
continuity  with  the  arc  corresponding  to  x  <  0.  Once  this  continuity 
is  perceived,  it  is  immediately  rationalized:  the  plane  (y,y)  cuts 
the  torus,  leaving  part  of  it  visible  above  it,  while  the  rest  remains 
hidden  under  the  plane  on  which  the  section  is  projected. 

In  one  respect,  Fig.  3  may  seem  odd:  the  invariant  sections  do  not 
contain  the  equilibrium  L^.  But  this  peculiar  feature  is  just  another 
instance  of  the  fact  that  we  are  dealing  here  with  the  projection  of  a 
surface  onto  a  plane.  Had  we  made  the  projection  or.  the  plane  (x,y), 
we  would  have  found  that  section  curves  in  the  x-hyperplane  do  indeed 
go  around  the  equilibrium  (see  Fig.  A). 

Most  analytical  treatments  of  an  equilibrium  treated  so  far  were 
concerned  with  Hamiltonians  without  gyroscopic  terms  (linear  in  the 
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Fig.  A.  -  Projection  on  the  plane  (x,y)  of  the  section  in  the  energy 
manifold  C  =  3  by  the  x-hyperplane  going  through  L^. 
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momenta)  wherein  the  equilibrium  constitutes  a  minimum  on  the  potential 
surface.  But,  in  the  restricted  problem,  and  are  maxima  on 

the  potential  surface,  whose  elliptic  character  in  stability  is  conferred 
by  the  Coriolis  forces  due  to  the  rotation  of  the  synodical  frame  of 
reference.  One  has  to  expect  that  the  global  structures  of  the  phase 
space  differ  substantially  on  the  energy  manifolds  according  as  to  whether 
they  lie  above  (C  <  3)  or  below  (C  >  3)  the  equilibrium  configur¬ 
ation.  We  illustrate  this  change  of  patterns  in  the  normalized 
system  cV  in  Fig.  5  where  we  present  the  projection  of  the  section 
of  the  energy  manifold  C  =  2.9  by  the  x-hyperplane  passing  through 
L,  . 

Since  we  have  in  analytical  form  the  general  solution  for  the 
principal  part  cMy  we  can  without  much  labor  observe  the  evolution  of 
invariant  sections  with  changes  in  the  Jacobi  constant.  By  way  of 

9 

illustration,  we  followed  the  projection  on  the  plane  (y ,y)  of  the 
section  by  the  x-plane  going  through  of  the  energy  manifold  for 

tue  Jacobi  constant  going  from  C  =  3.0001015  to  C  =  2.9.  On  each 
manifold  we  concentrated  on  the  ordinary  family  of  periodic  orbits 
0(25/2)  corresponding  to  the  rotation  number  p  =  25/2.  The  sequence 
of  curves  (see  Fig.  6)  begin  with  two  points:  these  are  the  intersec- 

n 

Lion  points  of  the  orbit  in  the  natural  family  Q  = £,  out  of  which 

1  A 

branch  the  ordinary  family.  Then  for  a  while,  the  invariant  section 
is  made  up  of  two  disjoint  ovals;  such  is  the  case  at  C  =  3.  Later 
on,  the  components  touch  one  another,  and  thereon  the  section  takes 


Fig.  5.  -  Projection  on  the  plane  (y,y)  of  the  section  in  the  energy 
manifold  C  =  2.9  by  the  x-hyperplane  going  through  L^.  (Normalized 
system) . 
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the  shape  of  two  hinged  crescents.  There  will  come  a  value  of  C  at 
which  the  crescents  will  separate;  each  closed  component  will  contract 
and  eventually,  at  C  =  2.559,  they  will  collapse  each  into  a  point 

g 

corresponding  to  the  orbit  in  the  natural  family  0  ^  ^  on  which  the 

ordinary  family  terminates. 

The  same  pattern  is  found  for  the  sections  in  the  ordinary  family 
of  periodic  orbits  <9(12/1)  (see  Fig.  7).  We  know  that  the  family 
originates  at  C  =  2.960750380  from  an  orbit  in  the  natural  family 

s 

£ and  that  it  terminates  at  C  =  2.635008320  on  an  orbit  in  the  same 
family.  By  the  way,  in  the  full  restricted  problem,  the  ordinary  two- 
parameter  family  of  periodic  orbits  associated  with  the  commensurability 
12/1  dissolves  into  a  natural  family  which  is  a  two  lane  bridge  from 
£S^  to  £ Of  the  curves  drawn  in  Fig.  7,  there  remain  only  finite 
sets  of  isolated  fixed  points,  half  of  them  being  elliptic  (for  the 
stable  orbits  in  the  bridge)  and  the  other  half  being  hyperbolic  (for 
the  unstable  orbits).  In  this  respect,  Fig.  7  which  refers  to  the 
truncated  normalized  system  cH  should  be  compared  with  the  corresponding 
figure  published  elsewhere  (Deprit  and  Rabe  1969)  for  the  full  restricted 
problem  of  three  bodies. 

So  far  we  have  considered  invariant  sections  in  the  map  defined 
by  the  original  Cartesian  coordinates.  It  is  not  without  interest  to 
analyze  them  in  the  normalizing  coordinates  ( ,1*  ,1*)  themselves. 
Thus  in  Fig.  8,  we  have  plotted  for  three  values  of  the  Jacobi  constant 
C  the  values  of  1*  and  l?)  of  the  points  at  which  the  periodic 
orbits  in  the  ordinary  family  0(25/2)  intersect  the  x-hypcrplane 
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Fig.  b.  -  Evolution  of  the  projection  on  the  plane  (y,v)  of  the 
sections  in  the  ordinary  family  0(25/2). 


I'ig.  7.  -  Involution  of  the  projection  on  the  plane  (y,y) 
sections  in  the  ordinary  family  of  periodic  orbits  f?(12/l) 
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going  through  L^.  Actually  Fig.  8  is  but  the  transform  of  Fig.  6 

by  the  normalizing  mapping  (x.y.p^p^)  ->  (4>* ,4>*,I*, I*) .  At 

C  =  3.0001015,  i.e.  the  energy  level  at  which  the  ordinary  family 

£ 

bifurcates  from  the  natural  family  £^,  the  section  would  consist  of 
two  lines  =  constant:  the  family  0( 25/2)  reduces  there  to  an 
orbit  of  long  period.  Somewhere  between  C  =  2.995  and  C  =  2.9, 
the  two  branches  of  the  intersection  curve  made  a  contact,  and  there¬ 
after  the  curve  reduces  to  a  closed  oval.  As  it  expands  in  the  direc¬ 
tion  of  the  angles  if*,  it  will  eventually  touch  its  image  in  the 
congruence  modulo  2tt  which  defines  the  angle  <f*.  From  this  level 
of  energy  onward,  the  intersection  curve  will  split  into  two  separate 
arcs.  The  evolution  terminates  when  these  arcs  have  become  two  straight 
lines  4* jjf  =  constant.  We  are  then  at  the  energy  level  C  =  2.559  at 
which  the  ordinary  family  0(25/2)  collapses  into  a  short  period  orbit 
circuited  25  times. 

The  preceeding  explorations  cause  us  to  question  an  empirical  defini¬ 
tion  given  lately  for  the  rotation  number  p.  Given  a  dynamical  system, 
usually  not  integrable,  the  section  of  an  energy  manifold  by  a  hyper¬ 
plane  is  determined  by  numerical  integration:  on  an  orbit  starting 
from  a  selection  of  initial  conditions,  the  program  marks  the  points  at 
which  the  hyperplane  is  penetrated.  In  the  absence  of  any  qualitative 
information  about  the  structure  of  the  phase  space  around  this  orbit, 

(1967)  proposes  to  define  the  rotation  number  as  the  average 
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Fig.  8.  -  Evolution  of  the  projection  on  the  torus  (,;*,;*)  of  the 
sections  in  the  ordinary  family  C?( 25  /2)  by  the  x-hvperpl ane  going 
through  L,. 
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angle  between  successive  pointers  to  the  intersection  points  from 
what  seems  to  be  the  fixed  point  inside  the  island  circumscribed  by 
the  "invariant  curve".  Let  us  explain  now  on  the  principal  part 
cV  for  the  restricted  problem  at  how  this  empirical  definition 

fails  to  be  satisfactory.  At  the  energy  level  C  =  2.560  a  bit  after 
the  family  branched  out  of  a  short  period  orbit  traveled  25  times, 
the  orbits  of  0(25/2)  have  25  well  developed  loops,  and  they  inter¬ 
sect  the  section  plane  in  25  points.  But  it  takes  two  full  rotations 
around  the  fixed  point  which  is  the  trace  of  the  short  period  orbit 
before  the  crossings  come  back  on  the  first  intersection  point.  Hence, 
according  to  Contopoulos ,  the  rotation  number  would  lie  in  the  vicinity 
of  25/2.  But,  as  C  increases,  the  loops  contract,  so  that  fewer  of 
them  are  crossing  the  section  hyperplane.  For  instance  at  the  level 
C  =  2.9,  we  count  only  16  intersections  on  the  orbits  of  0( 25/2)  after 
two  full  rotations:  according  to  Contopoulos*  definition,  the  rotation 
number  would  have  decreased  now  to  16/2.  Eventually  the  family  approaches 
a  long  period  orbit  traveled  twice.  A  little  before  C  =  3.0001015  where 
the  bifurcation  takes  place,  the  plane  section  presents  only  two  crossing 
points,  in  which  case  Contopoulos*  rotation  number  cannot  be  unambiguous¬ 
ly  estimated. 
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3.  The  resolution  power  of  a  normalization 

How  does  the  portrait  in  phase  space  for  the  normalized  system 
oV  compare  with  that  of  the  full  system?  In  order  to  answer  that 
question,  it  might  seem  logical  to  begin  by  inquiring  about  the  con¬ 
vergence  of  Birkhoff's  normalizing  transformations.  The  few  precise 
results  mathematicians  have  so  far  established  look  discouraging. 

We  know  that,  in  an  admissible  set  of  power  series  provided  with  a 
natural  topology,  the  subset  of  Hamiltonians  for  which  Birkhoff's  normali¬ 
zations  are  divergent  is  everywhere  dense  (Siegel  1941)  ,  whereas  the 
subset  of  Hamiltonians  for  which  these  transformations  are  convergent 
is  of  first  category  in  the  sense  of  Baire  (Siegel  1954).  However, 

Siegel's  topology  can  be  weakened  sufficiently  without  yet  losing 
a  certain  tightness  so  that  the  set  of  Hamiltonian  systems  for  which 
Birkhoff's  normalization  converges  appears  as  everywhere  dense  (Contopoulos 
1963)  . 

Such  theorems  show  that  so  far  the  mathematical  problem  of  con¬ 
vergence  is  not  yet  of  physical  significance.  From  a  purely  experimen¬ 
tal  point  of  view,  we  have  learned  to  use  normalization  as  a  tool,  much 
as  a  physicist  uses  a  microscope.  The  ultimate  structure  in  the  phase 
space  around  the  equilibrium  is  yet  beyond  our  reach.  Numerical 
integration  of  the  original  equations  of  motion  will  provide  a  detail 
here  and  there,  but  it  will  leave  us  totally  ignorant  of  some  global 
structures.  For  instance  we  might  be  interested  in  locating  some 
invariant  tori  of  quasi-periodic  motions  surrounding  the  equilibrium. 
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Let  us  assume  that  we  are  in  the  conditions  which  guarantee  their 
existence  for  the  full  Hamiltonian  In  this  case,  we  could  act 

under  the  assumption  that,  if  the  transformation  is  carried  to  a 
sufficiently  high  degree,  the  normalized  system  cK  will  possess  tori 
in  the  domain  of  phase  space  that  we  analyze,  and  that  the  perturbing 
residual  jP  will  not  wipe  them  all  from  the  phase  space  of  the  full 
system. 

Normalization  will  reveal  some  of  the  features  in  the  phase 
space  around  the  equilibrium.  Of  course  on  the  whole  it  will  simplify 
them.  Where  it  tells  the  story  of  two  basic  frequencies  so  well  locked 
in  a  resonance  that  they  generate  an  ordinary  family  of  periodic  orbits, 
we  have  to  expect  that,  in  the  full  problem,  the  ordinary  family  breaks 
down  into  isolated  natural  families  of  periodic  orbits,  some  of  them 
being  of  the  elliptic  type  and  thus  surrounded  by  "islands"  indicative  of 
tori  of  quasi-periodic  orbits,  others  being  of  the  hyperbolic  type  and 
acting  as  epicentfes  in  the  middle  of  instability  zones.  But  this  dis¬ 
tribution  of  stability  islands  and  instability  faults  is  very  likely 
inscribed  in  annulus  bordered  by  fringes  of  quasi-periodic  notions 
with  nonzero  measure.  The  normalization  is  unable  to  resolve  these 
fine  details;  yet  through  the  blurring,  an  experienced  investigator 
might  gain  a  feeling  for  what  complexity  there  is  in  the  full  problem. 

A  normalization  truncated  after  the  first  few  orders  may  be  out¬ 
rageously  simplifying.  The  higher  the  order  of  truncation  is,  the  more 
significant  the  approximation  could  be.  Compare  for  instance  Fig.  12, 
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13,  and  1  of  Gustavson  (1966)  where  a  progress  for  order  6  to  8 
results  in  bringing  forth  more  patently  new  fixed  points  in  the 
sections,  i.e.  natural  families  in  the  phase  space.  Of  course,  like 
with  any  asymptotic  series,  there  comes  an  optimum  beyond  which  the 
resolution  power  of  the  normalization  loses  its  significance.  And 
it  is  likely  that,  even  before  the  optimum  in  truncation  is  reached, 
numerical  accuracy  will  have  failed.  For  the  normalization  will  have 
been  obtained  by  so  many  operations  on  the  intermediate  results  that 
the  high  order  coefficients  in  the  power  series  will  have  lost  most, 
if  not  all,  of  their  relative  accuracy;  besides  coordinates  and  velocities 
obtained  by  adding  a  large  number  of  small  terms  soon  cease  to  be 
meaningful. 

Whereas  the  mathematician  spontaneously  looks  for  bounded  domains 
in  which  trajectories  can  be  guaranteed  to  flow  for  ever  without  leak¬ 
ing  or  for  uniform  vicinities  in  which  orbits  are  predicted  to  stay 
for  ever  close  to  one  another,  astronomers  and  physicists  accommodate 
themselves  with  more  realistic  and  thus  less  stringent  requirements. 

In  order  to  build  ephemerides  for  Trojan  planets  from  an  analytical 
theory,  it  is  sufficient  to  be  certain  that  an  intermediate  orbit 
provided  by  a  normalizing  transformation  about  will  stay  for  a 

sufficiently  long  time,  although  not  for  ever,  reasonably  close  to 
an  approximately  quasi-pe riodic  motion  of  the  restricted  problem  of 
three  bodies.  In  this  respect  we  see  two  reasons  why  the  normalized 


model  mav  fail: 
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(i)  As  we  carry  the  normalization  up  to  an  assigned  order, 
truncation  errors  will  steadily  decrease  the  accuracy  of  the  normali¬ 
zation  as  the  normalized  action  integrals  move  away  from  their  singu¬ 
lar  values  at  the  equilibrium. 

(ii)  But  also  as  the  order  of  the  normalization  goes  to  infinity, 
we  do  not  expect  that  the  normalizing  transformations  converge  uni¬ 
formly  in  a  neighborhood  of  the  equilibrium.  Rather  we  shall  now  show 
that  there  are  regions  where  the  structures  of  the  phase  space  for  the 
full  problem  differ  significantly  from  the  ones  predicted  by  the  normali¬ 
zation.  In  these  regions  we  shall  speak  of  model  errors. 

As  we  plan  to  test  the  normalizing  approximation  where  we  expect 
it  to  be  at  its  worst,  we  focused  our  attention  on  the  ordinary  families 
of  periodic  orbits  revealed  by  the  normalization.  Another  test  would 
have  been  to  compare  for  the  same  initial  conditions  the  solution  of 
the  restricted  problem  and  the  quasi-periodic  motion  of  its  normalization. 

Actually  the  test  on  ordinary  families  of  periodic  orbits  has  been 
done  under  very  convincing  conditions.  Four  years  before  we  even  con¬ 
templated  normalizing  the  restricted  problem  at  ,  we  began  computing 
invariant  sections  in  the  neighborhood.  The  differential  and  variational 
equations  were  integrated  simultaneously  by  recurrent  power  series 
(Deprit  and  Price  1965)  so  as  to  assure  at  least  twelve  significant 
figures  for  the  orbit  and  nine  significant  figures  for  the  variations 
for  at  least  several  hundred  units  of  time.  At  the  level  C  =  3.00007, 
we  obtained  at  first  several  curves  strongly  suggestive  of  quasi-periodic 
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motions;  then  we  began  looking  for  stable  periodic  orbits  belonging 
to  natural  families  and  thus  causing  swirls  to  appear  in  the  general 
stream  of  quasi-periodic  motions.  To  our  dismay,  all  we  could  get  was 
several  sets  of  periodic  orbits  with  equal  periods  and  the  stability 
index  Tr(T)  equal  to  2  up  to  eight  decimal  tigures.  Moreover,  in 
each  case,  the  normal  variation  n£  solution  of  Hill's  homogeneous 
equation  for  the  initial  conditions  n2(0)  =  0,  n2(0)  *  1,  and  the 

normal  variation  n^  solution  of  Hill's  nonhomogeneous  equation  for 
the  initial  conditions  n^(0)  *  n^CO)  =  0  checked  at  the  end  of  the 
period  the  relation 

n2(T)n3(T)  -  [^(T)  -  l]n3(T)  -  0, 

as  predicted  for  an  ordinary  periodic  orbit  (Deprit  and  Henrard  1968). 

Initial  conditions  for  these  remarkable  orbits  constitute  all  the 
entries  of  Table  I  but  the  last  two  ones.  The  integrations  have  been 
carried  out  at  two  different  seasons,  hence  unfortunately  a  slight  varia¬ 
tion  of  the  mass  ratio  which  reflects  in  the  entries  of  the  first  row. 

It  has  induced  small  alterations  in  the  corresponding  periods. 

The  abundance  of  seemingly  ordinary  families  of  periodic  orbits 
seemed  odd  to  ns  at  the  time  we  discovered  them.  For  in  the  opinion 
of  some  experts,  "ordinary  solut ions ...  are  never  found  after  extensive 
search  by  many  investigators"  (llray  and  Goudas  1967).  Normalization 
ot  the  restricted  problem  of  three  bodies  was  then  undertaken  to  eluci¬ 


date  these  results. 
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We  started  with  initial  conditions  which,  for  the  normalized 
model,  correspond  to  ordinary  families  of  periodic  orbits.  We  inte¬ 
grated  numerically  the  orbits  they  generate  in  the  restricted  problem 
and  stopped  the  calculation  when  we  reached  the  estimate  of  the  period 
predicted  by  the  normalization.  Then  either,  within  the  assigned  tolerance, 
the  final  state  was  equal  to  the  initial  one,  or  there  was  a  discrepancy. 

In  the  former  case,  the  orbit  was  held  to  be  periodic.  In  the  latter, 
by  trial  and  error,  the  initial  conditions  were  adjusted  so  as  to  pro¬ 
duce  a  periodic  orbit. 

If  its  characteristic  exponents  came  all  four  very  close  to  zero, 
then  we  restarted  the  search  for  periodic  orbits  along  the  invariant 
curve  of  the  ordinary  family  of  periodic  orbits  in  the  normalized 
problem.  This  was  a  cumbersome  and  time-consuming  investigation;  the 
usual  differential  corrections  to  converge  onto  a  periodic  orbit  fail  when 
the  characteristic  exponents  are  so  definitely  close  to  zero.  Thus 
we  had  to  satisfy  ourselves  with  a  few  orbits  for  each  invariant  case. 
Anyway,  first  we  recovered  in  this  manner  the  ordinary  families  that 
had  been  discovered  prior  to  the  normalization,  and,  then,  on  each 
of  them,  we  added  a  few  more.  It  is  reasonable  to  conclude  that,  in 
a  close  neighborhood  of  ,  within  the  accuracy  at  present  available 
by  careful  calculations  in  double  precision  the  structure  of  the  phase 
space  about  cannot  be  magnified  beyond  the  resolution  of  a  Birkhoff's 

normalization  of  the  thirteenth  order.  The  ordinary  families  of  the 
normalized  system  for  the  rotation  numbers  25/2,  35/3,  51/4,  64/5  and 
77/6  could  not  be  broken  down  into  natural  families  of  periodic  orbits. 
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Yet  further  away  from  the  equilibrium  in  the  phase  plane  (I-,!*), 
the  torus  of  the  normalized  system  which  is  filled  with  the  ordinary 
periodic  orbits  corresponding  to  the  rotation  number  13/1  is  washed 
away  from  the  restricted  problem.  In  its  place,  numerical  integration 
detected  long  and  thin  islands.  But  the  crossing  points  on  the  section 
curves  move  so  slowly  that  we  did  not  think  it  worthwhile  to  calculate 
the  complete  pattern.  In  one  case,  we  had  followed  a  trajectory  for 
1600  units  of  time  (i.e.  30000  solar  years)  ,  and  yet  we  had  not  even 
traveled  halfway  around  the  island.  Instead  we  applied  variational 
corrections  to  converge  onto  the  subsiding  natural  periodic  orbits. 

One  of  them  (see  the  last  two  lines  of  Table  I)  has  unstable  charac¬ 
teristic  exponents,  and  we  discovered  later  on  by  numerical  continua- 
tion  that  it  is  an  element  of  the  family  £  ^  on  Us  way  back  from  its 
maximum  Jacobi  constant.  The  other  natural  periodic  orbit  has  characteris¬ 
tic  exponents  of  the  stable  type;  we  also  established  later  on  that  it 

belongs  to  the  bridge  J(13S,14S)  of  periodic  orbits  connecting  an 

s  s 

element  of  £ ,  traveled  13  times  to  another  element  of  £ .  traveled 
4  4 

14  times. 

The  dissolution  of  the  family  £7(13/1)  could  be  rationalized 
on  two  grounds:  not  only  does  the  resonance  13/1  occur  in  the  normalized 
system  further  away  from  the  equilibrium  than  the  commensurabi lities 
previously  encountered,  but  also,  as  indicated  by  Contopoulos  (1967), 

the  width  of  an  island  at  the  resonance  m/n  is  proportional  to 

*s(m+n-4)  .  .  ,  , 

f  ,  where  e  is  an  estimate  of  the  distance  to  the  equilibrium. 

In  the  present  case,  the  sum  m  +  n  drops  to  its  minimum  value  when 
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the  rotation  number  passes  by  the  ratio  13/1,  and  therefore  in  an 
intuitive  way,  we  expect  significantly  wider  islands  along  the  in¬ 
variant  curve  of  0(13/1) . 

The  breakdown  of  an  ordinary  family  into  a  pair  of  natural 
families  with  characteristic  exponents  of  opposite  type  in  stability 
is  but  an  illustration  of  familiar  statements  of  Poincare.  What  it 
tells  also  is  that  the  normalization  is  only  a  local  tool,  and  that 
it  fails  to  reveal  the  global  structure  of  the  phase  space.  Indeed, 
in  the  normalized  approximation,  the  ordinary  family  0(13/1)  consti- 
tutes  a  two-parameter  bridge  connecting  a  natural  orbit  of  £  ^  and 
a  natural  orbit  of  X’f  travelel  thirteen  times.  But,  in  the  restricted 
problem,  the  two  natural  families  which  survive  the  destruction  of 
0(13/1)  have  certain  terminations  modified. 

Our  findings  at  Lhe  energy  level  C  =  3.00007  are  presented  in 
Fig.  9. 

A  similar  analysis  has  been  undertaken  at  the  energy  level  C  = 

2.90.  We  find  there  the  same  pattern.  Close  to  the  fixed  point  belong- 

g 

ing  to  the  short  period  orbit — an  element  of  the  natural  family  £  ^ — the 
ordinary  families  of  the  normalized  system  subsist  in  the  restricted 
problem.  Such  is  the  case,  for  instance,  of  the  family  0(35/3). 

But  at  a  greater  distance,  the  invariant  curves  of  the  ordinary  families 
break  down  into  interspersed  subsets  of  isolated  fixed  points  correspond¬ 
ing  to  natural  families  of  periodic  orbits.  The  invariant  curve  of 
the  ordinary  family  0(12/1)  gives  way  to  the  elliptic  fixed  points 
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(the  small  circles  on  Fig.  10)  of  a  stable  periodic  orbit  and  the 
hyperbolic  fixed  point  (the  crosses  on  Fig.  10)  of  an  unstable  periodic 
orbit;  they  belong  respectively  to  the  stable  and  the  unstable  lane 
of  a  bridge  3>(12S,12S)  connecting  a  short  period  orbit  at  C  = 
2.960750380  traveled  12  times  to  the  short  period  orbit  at  C  = 
2.635008320  also  traveled  12  times  (Deprit  and  Rabe  1969).  The 
ordinary  family  (7(25/2)  dissipates  likewise,  making  place  to  the 
stable  and  unstable  lanes  of  a  natural  bridge  JJ(25S,2L)  connecting 
the  short  period  orbit  at  C  =  2.559466612  traveled  25  times  to  the 
long  period  orbit  at  about  C  =  3.0003  traveled  twice.  Also  the  curve 
of  (7(13/1)  survives  only  in  the  elliptic  fixed  points  imprinted  by 

the  orbit  in  the  stable  leg  of  the  bridge  5(13S,14S),  and  the  hyper- 

0 

bolic  fixed  points  of  the  unstable  orbit  in  the  family  £;  on  its 
way  from  a  maximum  C  to  its  termination  on  the  short  period  orbit 
at  C  =  2.501429345.  Beyond  this  invariant  curve,  the  normalization 
truncated  beyond  order  13  is  no  longer  a  valid  approximation.  In  the 
shaded  area  ol  Fig.  10,  we  have  marked  by  x's  the  intersection  points 
ot  the  unstable  periodic  orbit  in  the  return  leg  of  tb<_-  bridge  J>(18S,19S), 
although  we  have  not  tried  to  elucidate  the  structure  of  the  phase  space 
in  that  remote  part  of  the  energy  manifold. 

In  this  respect  we  cannot  leave  our  readers  unaware  of  the  over¬ 
abundant  amount  of  numerical  integrations  required  to  produce  but  the 
handful  of  conclusions  we  came  to  so  far.  As  an  illustration  we  have 
blown  up  in  Fig.  11  a  minute  area  of  the  cross  section  plotted  in 
Fig.  10.  This  is  a  rectangle  covering  the  dissipated  invariant  curve 
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relative  to  0(25/2)  between  two  successive  hyperbolic  fixed  points 
belonging  to  the  two-lane  bridge  2?(25S,2L).  The  closed  curve  around 
the  elliptic  fixed  point  demanded  that  the  numerical  integration  be 
extended  over  5723  units  of  time  (i.e.  about  10,000  solar  years),  yet 
keeping  fixed  10  decimal  figures  of  the  Jacobi  constant. 

At  the  level  C  =  3.00007,  the  ordinary  family  0(25/2)  of  the 
normalized  system  seemingly  subsists  in  the  restricted  problem,  at 
least  within  the  constraints  of  an  accuracy  limited  to  the  double 
precision  arithmetic  of  an  IBM  7094;  but  at  C  =  2.9,  it  has  dissolved. 
The  dissolution  is  progressive  as  the  Jacobi  constant  recedes  from  its 
value  C  =  3  at  the  equilibrium  state.  As  the  natural  family  J(25S,2L) 

g 

moves  from  its  bifurcation  out  of  £  we  noticed  that  the  stability 
index  at  first  grows  larger  than  2  along  the  unstable  lane  of  the  bridge, 
then  decreases  toward  2  and,  starting  with  C  =  2.99,  it  stays  equal 
to  2  up  to  nine  decimal  figures  (see  Fig.  12).  At  this  point  the  natural 
t.imily  has  entered  the  zone  around  where  the  restricted  problem 

cannot  be  detached  from  its  normalizing  approximation. 

The  ne ighborhoods  of  the  equilibrium  that  we  considered  often 
times  in  this  analysis  are  of  course  defined  in  the  plane  of  normalized 
actions  (I*,  I*).  Mapped  back  into  the  original  Cartesian  phase  space, 
the  orbits  we  had  under  study  usually  depart  significantly  from  the 
textbook  concept  of  an  infinitesimal  neighborhood  of  I...  Bv  wav  of 
illustrating  this  contrast,  we  have  drawn  in  Fig.  13  the  unstable  orbit 
of  the  bridge  J>(12S,12S)  that  we  discovered  on  the  energy  manifold 
C  = 


2.9 


Fi^.  12.  -  Stability  curve  for  the  stable  rnd  -instable  lanes  of  the 
bridge  J»(  2  5S  ,21.)  showing  the  progressive  blu:  ring  from  a  natural 
family  into  the  ordinary  family  0(25/2). 
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Conclusions 

The  relations  between  a  dynamical  system  having  an  equilibrium 
and  its  truncated  Birkhoff's  normalization  around  that  critical  element 
may  stir  as  much  controversy  as  the  existence  of  the  so-called  third 
integral  in  galactic  dynamics.  As  a  matter  of  fact  the  two  formalisms 
bear  on  the  same  technique  which  is  to  approximate,  of  course  in  a 
conveniently  unspecified  topology,  a  dynamical  system — integrable  or 
not — by  one  whose  variables  are  obviously  separable.  The  technique 
of  the  third  integral  consists  in  building  formally  an  integral  that 
is  both  independent  of  the  energy  and  in  involution  with  it;  from 
thereon  one  relies  on  Liouville's  theorem  to  provide  the  capability  of 
constructing  a  system  of  separable  coordinates.  On  the  other  hand, 
Birkhoff's  normalization  aims  straightforwardly  at  constructing  these 
separable  coordinates,  in  such  a  way  that,  the  angle  coordinates  being 
ignorable,  the  action  momenta  turn  out  to  be  the  independent  integrals 
in  involution  sought  by  the  technique  of  the  third  integral. 

From  the  mathematical  standpoint,  the  basic  concern  so  far  has 
been  with  deciding  whether  or  not  the  third  integral  exists,  or  equivalently 
whether  Birkhoff’s  normalizing  transformations  form  a  convergent  sequence. 
Whatever  the  answer  to  that  question  may  be,  and  it  may  well  be  in  the 
negative,  one  cannot  but  wonder  at  the  many  and  deep  physical  insights 
these  mathematically  crippled  techniques  provide  in  astronomy.  In  a 
way  numerical  experimentations  in  celestial  mechanics  compel  us  to 
adopt  the  kind  of  pragmatism  that  pervades  observational  astronomy. 
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Birkhoff's  normalization  acts  as  a  lens  focused  on  a  dynamical  system. 

It  helps  in  penetrating  the  complexity  beyond  the  first  order  appear¬ 
ances.  Yet  its  resolution  power  is  limited.  There  comes  a  threshold 
beyond  which  the  present  computer  equipment  does  not  allow  us  to  separate 
the  normalized  image  from  the  real  phase  space.  Fortunately,  at  least 
in  the  case  we  have  had  under  special  investigation,  to  learn  that  there 
is  a  hard  core  in  which  double  precision  arithmetic  cannot  divide  the 
restricted  problem  from  its  normalization  truncated  beyond  order  13  is 
not  unwelcome  news.  At  least  it  offers  a  way  of  substantially  improving 
the  theory  of  Trojan  planets  by  building  it  on  Birkhoff's  normalization. 

On  the  whole  the  present  study  v.  .  ents  numerically  the  views  held 
by  Contopoulos  and  others  (see  for  instance  Walker  and  Ford  1969).  For 
librations  around  the  equilibrium  with  sufficient^  small  amplitudes,  >  t 
orbits  are  not  influenced  by  resonances  in  a  way  that  could  presently  be 
detected  by  accurate  integrations  carried  over  a  long  time;  yet  one  is 
already  able  to  detect  exceptional  motions  influenced  by  isolated  resonances. 
Mathematically  speaking  these  motions  are  densely  populating  the  area 
covered  by  the  well  behaved  majority;  however  Birkhoff's  normalization 
filters  only  a  few  of  the  exceptional  motions,  namely  those  corresponding 
to  the  smallest  resonances,  and  thus  isolates  them.  We  have  shown  here 
that  numerical  integration,  however  accurate,  could  do  no  better  than 
normalization  in  this  respect. 
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Appendix  1 

The  program 

The  basic  information  on  how  we  represent  Poisson  series  (3-3) 
in  the  computer  has  been  given  elsewhere  (Deprit  and  Rom  1968  a,  1968  b) . 
We  concentrate  here  on  the  specific  program  which  implements  Birkhoff's 
normalization. 

'file  function  BCOEF(N.M)  computes  in  double  precision  arithmetic 
the  binomial  coefficient  (!^)*  Dr.  David  Walkup,  our  colleague  at 
B.S.R.L.,  is  the  author  of  this  program. 

Before  we  describe  the  next  subroutines,  let  us  agree  on  some 
conventions.  We  operate  exclusively  with  Poisson  series  of  the  type 
(3-3);  elements  of  a  series  are  of  the  form 

n  ^  ii/^ 

cosCm^  +  m0;0  +  m^), 

n |  n9 

x^  x0“x3‘  sin(nij;^  +  m9:0  +  m^;^). 

I  he  first  polynomial  argument  is  not  used,  the  second  stands  for  /T  , 
and  the  third  for  >  S  ;  the  first  trigonometric  argument  refers  to  £, 
the  second  refers  to  s,  and  the  third  is  not  used. 

The  subroutine  PDS(A,B,J)  differentiates  the  Poisson  series  A 
with  respect  to  the  polynomial  argument  Xj  and  stores  the  result  as 
the  Poisson  series  B. 

The  subroutine  PDJS (A,B  ,  . ,y)  applies  to  the  Poisson  series  A 


tiie  differential  operator 
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a 


0  t —  + 

3  c  ^ 

2 


3 


and  stores  the  result  as  the  Poisson  series  B. 

The  subroutine  PR0D3(A,B ,C, a)  stores  as  the  Poisson  series  C 
the  result  of  the  operation  C  +  aAB. 

These  three  commands  enable  us  to  compute  the  Jacobian  determinants 

3 A  3B _ _3A_  _3B_  3A_  _3B_  3A  3B 

3$^  3x2  3x2  3<J>^  ’  34>2  3x^  3x^  3<f>2 


But,  when  we  compute  the  Poisson  bracket  of  two  functions  in  L,  S,  1 
and  s,  we  should  not  overlook  that  x2  stands  for  /L  and  not  for 
L,  and  that  x^  stands  for  /"s’  and  not  for  S.  Therefore 

3 ( A ; B )  =  1  1_  3 (A,B) 

3 (£,])  2  x2  3(^1,x2)  ’ 

3 (A,B)  m  11_  3 (A,B) 

3(s,S)  2  x^  3(^2  »xj) 

Therefore  we  need  the  subroutine  DVLS(A,J)  which  replaces  the  Poisson 
series  A  by  the  series  A/xj. 


The  various  partial  derivatives  of  a  d'Alembert  series  and  their 
products  do  not  necessarily  possess  the  d'Alembert  characteristic.  But 
the  exceptional  terms  in  the  products  —  and  —  --g  cancel  one 
another  in  the  formation  of  the  Jacobian  determinant  •  1°  the 

computer,  should  we  proceed  in  integer  arithmetic,  all  the  exceptional 
terms  would  disappear  in  the  ’acobian  which  thus  would  possess  the 
d'Alembert  characteristic  exactly.  But,  as  we  proceed  only  in  double 


♦ 
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precision,  the  expected  cancellations  are  missed  by  relatively  small 
amounts.  The  subroutine  DALPIC(A,B)  remedies  this  deficiency:  it 
scans  the  elements  of  A  checking  for  each  one  such  that  |m^| 

|m9|  1  n3,  n^  -  (mod.  2)  and  (mod.  2).  The  terms  which 

satisfy  these  properties  are  stored  as  the  Poisson  series  B. 

The  preceding  subprograms  are  called  in  the  subroutine 
DPBRAC( /A/ , /B / , /C/) .  The  Poisson  series  A  and  B  are  assumed  to 
have  the  d'Alembert  characteristic,  the  resulting  Poisson  series  C 
is  the  d'Alembert  series 

r  m  '(A,B)  3(A,B) 

C  3(e,L)  Hs,S)  • 

i.e.,  the  Poisson  bracket  (A;B),  of  the  functions  A  A(£,s,L,S) 
and  B  B(f,s,L,S)  in  the  phase  space  (£,s,L,S). 

DP BRAC  is  the  main  ingredient  oi  the  subroutine  LTRAN(F,W ,P ,Q ,N , /ANS/)  . 
It  is  assumed  that  all  Poisson  series  entering  the  triangle 


(1) 


.(1) 

'1 


f 


(2) 

0 


(1)  (2)  (3) 

2  1  0 


liave  been  defined  up  to  and  including  those  belonging  to  the  (n-l)-th 
row.  Note  that  these  elements  are  labeled  in  FORTRAN  as  follows 
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F(l,l) 


F(l,2)  F(2 ,1) 

F(1 ,3)  F(2 ,2)  F(3, 1) 

F(l,4)  F(2,3)  F(3 ,2)  F(<  ,1) 


Also  it  is  assumed  that  the  Poisson  series  W,  ,W„,...,W  have  received 

12  n 

a  FORTRAN  name,  W(l)  ,W(2)  , . . . ,W(N) .  The  subroutine  computes  the  element 
f^  from  the  rule 

q 


(p),  V  ,Qyt(p-i,.u  \  +  f(p-i> 

Q  ’Vi'  Vi 


and  stores  it  as  the  Poisson  series  ANS.  One  should  not  overlook  that 
LTRAN  assumes  that  input  and  output  series  have  the  d'Alembert  charac¬ 
teristic;  it  cannot  be  regarded  as  a  general  subroutine  to  compute  any 
element  in  the  triangle  of  a  Lie  transform  under  any  canonical  trans¬ 
formation. 


BIRKOF(NMAX)  essentially  performs  a  cycle  to  generate  from  the 

(1)  (NMAX) 

sequence  J{Q  , .  .  .  ^NMAX  the  sequence  , . . .  ,  that 

is  the  normalized  Hamiltonian,  and  the  sequence  W,  ,W_  , .  . .  ,WKTU  that 
defines  the  normalizing  transformation.  The  N-th  cycle  begins  by  set¬ 
ting  equal  to  zero;  it  is  accomplished  by  declaring  that  — whose 

FORTRAN  designation  is  W(N+1) — is  the  Poisson  series  in  formation  at 
the  bottom  of  the  stack.  As  long  as  no  element  is  entered  in  ,  it 

will  represent  the  null  d'Alembert  series.  Then  the  provisional  ele- 
~(1)  ~(2)  ~(N) 

ments  » •  •  • »  are  formed  by  calling  recursively  LTRAN; 

~(N) 

they  are  stored  in  the  disk,  except 


which  is  kept  in  core  in 
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the  temporary  series  H.  The  subroutine  COLECT  transcribes  in 

(N)  ~"(N) 

the  nontrigonometric  terms  of  ,  as  indicated  by  the  value 

(0,0, 0,0)  given  initially  to  its  argument  TPART  (N.B.  Nontrigono¬ 
metric  terms  in  Poisson  series  of  type  (3,3)  are  characterized  as 
having  the  last  four  flags  of  their  label  word  equal  to  zero).  At 
this  point  we  have  obtained  a  new  component  of  the  normalized  Hamil¬ 
tonian. 

The  computation  of  begins  by  storing  in  P  the  purely  trigo¬ 

nometric  part  »£N)  by  means  of  the  subroutine  SCREEN.  Then  IPDJS 
accomplishes  the  quadrature 

W  =  — i -  P 

N  D  (  X  ,  c ) 

and  stores  the  resulting  Poisson  series  in  W^.  It  is  punched  in 

hexadecimal  (HEXOUT)  and  it  is  written  on  the  disk  (disk  address 

stored  in  W(N+1)).  Unless  the  program  has  reached  the  last  component 

WN  to  be  obtained,  the  triangle  elements  , . . .  ,  are 

then  completed  by  subtracting  from  each  of  them  the  periodic  part  P 

(N  ) 

of  3f0  .  The  N-step  of  the  cycle  ends  on  a  reading  of  the  clock  to 

show  the  time  spent  in  the  operations  at  that  stage. 

Hie  subroutine  BIRKOF(NMAX)  is  escorted  by  another  program  which 
produces  the  Lie  transform  of  the  coordinates  x,y  and  the  velocities 
x,y.  Obviously  it  is  built  around  a  cycle  calling  recursively  the  sub¬ 


routine  LTRAN. 


Appendix  2 


Listings 
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FORTRAN 

0001 

0002 

0003 

0004 

0005 

0006 

00C7 

oooe 

0009 
0010 
0011 
0012 
0013 
0014 
0015 
CO  1 6 
CO  1  7 
CO  1 8 
0019 
0C20 
0021 
C022 
0023 
002* 
0025 
C  026 
C027 
0028 
0029 
CC  30 
CC  3  1 
0032 
0033 
CO  3* 

CC  35 
0036 
OC  3  7 
CO  38 
C039 
C0*0 
CO*t 
C  0*2 
C043 
G044 
C045 
CC-fc 
C04  1 


IV  MODEL  44  PS  VERSION  2*  LEVEL  1  OATE 

SUBROUTINE  B I RKOF ( NMAX 1 
REAL'8  LAMROA, SIGMA 

COMCN/H  AM/LAMBDA,  SIGMA,  Ft  20, 20  )»W<  20) 
INTEGER  P,C.TPARTI4)/4'0/ 

CC  ICO  S* l , NM A X 
CALL  DEFINEtWlNM)  ) 

CALL  GETIMEtll) 

C«N-1 

CC  10  J«1,N 

CALL  LTRAN|F,W, J,Q,?0,H) 

IF(J.EC.N)  GC  TO  12 
CALL  SEEK(C) 

CALL  ViRITEPIH.FI  J*1,CM), -1,0) 

CALL  ERASE(H) 

10  C«C- l 

12  CALL  CCL ECT I H , F I N* 1 »  1  )  » TPART ) 

CALL  LlST3IF(N*l,l)»*H,,-N) 

CALL  HE XCUT (FIN' 1*1)  ) 

CALL  SCREEMH.P,  TPART) 

CALL  ERASEIH) 

CALL  IPDJSIP,WN,LAMBDA,SIGMA) 

CALL  HE  XCUT I NN ) 

CALL  SEEK(C) 

CALL  ERASE ( W I N ♦  1  )  ) 

CALL  WRITEPIwN.WlNM  ),-l,200*N) 

CALL  ERASE1WN) 

IFIN.EC.l)  GC  TO  101 
IF  IN. EC. NMAX)  GO  TO  101 
Nl  *N- l 
C  =  N 

CC  20  J*t,Nl 
J 1 « J>  1 

CALL  CEFINEITEMP) 

#  CALL  ACCM1 ( P, TEMP,-1 .00) 

CALL  ACUM1IHC,  Jl  )  ,TEMP,1.0D0) 

CALL  SEEKfF(C,Jl)  ) 

CALL  mRIIEPI TEMP,F(g,Jl),-l,0) 

CALL  ERASEITEMP) 

2C  C*C-l 

101  CALI  GETIMEI  12) 

ETWI2-I  1 

ETM«(ET*/1S2C0. )/60. 

PRINT  40,N,ETM 

40  F0RMATI5X,»TIME  FOR', 13,'  ■  *,F7.2> 

100  CALL  ERASEIP) 

RETURN 

ENT 
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FORTRAN  IV  NCDEL  44  PS  VERSION  2,  LEVEL  l  OATE 


OOOl 

SUBROUTINE  0P6RACI /F/ , /G/ , /PB/ ) 

0002 

REAL*8  LS(2,2>/l.D0,0.,0.,i.D0/ 

0003 

CALL 

CE  F I NE ( P8 1 

0004 

cc 

10  I-  1,2 

0005 

CALL 

CEFINEIP1I 

0006 

CALL 

P0JS(F,DF,LS(1,I)*LS(2, I), 

0007 

CALL 

PCS ( G , CG, I  +  1 ) 

0008 

CALL 

PRCC3 ( OF , DG,  P 1 ,  •  500 ) 

0009 

CALL 

ERASE ( OG) 

0010 

CALL 

ERASE  I  OF ) 

con 

CALL 

POS  IF,  of, mi 

0012 

CALL 

PCJS(G,CG,LS(l,n,LS(2,I), 

0013 

CALL 

PR003(0F,DGf P1,-.5D0) 

CO  14 

CALL 

ERASE ( OF ) 

0015 

CALL 

ERASEICG) 

0016 

CALL 

CRCERIPl) 

0017 

CALL 

0VLS(P1,M 

coie 

CALL 

ACUM1 (P1,PB, l.DO) 

0019 

10  CALL 

ERASE  1  Pi ) 

0020 

CALL 

CALPIC(PB.SPB) 

0021 

CALL 

ERASE ( P8) 

0022 

CALL 

PXCMPB.SPB) 

0023 

RETURN 

0024 

ENC 

FORTRAN  IV  PCOEL  44  PS  VERSION  2 ,  LEVEL  I  OATE 


0001 

SUBROUTINE  LTRAN(F,W,P,Q,N,/ANS/1 

0002 

REAL«8  BCOEF 

0003 

INTEGER  P,C 

0004 

0  1  PENS  I CN  FIN, N), WIN) 

0005 

Nj-en 

0C06 

CALL  CEFINEIANS) 

0007 

00  10  J 1 ■ 1 ,NJ 

0008 

CALL  0PBRAC(FIP,U-J1*2),WI jm),B 

0009 

CALL  ACUN1(B,ANS,BC0EFIC,J1-1) ) 

0010 

10 

CALL  ERASE  I B ) 

0011 

CALL  ACUN1IFIP, Q*2), ANS, 1.00) 

0012 

RETURN 

0013 

ENC 

FORTRAN  IV 

NCOEl 

44  PS  VERSION  2,  LEVEL  1 

CC01 

DOUBLE  PRECISION  FUNCTION  BCOEFIN, 

C002 

IBC-1 

0003 

IF (24P.LE.NIG0  TO  20 

CC04 

I  B  ■  N-  P 

C005 

GO  TC  30 

0006 

20 

I  8*P 

00C7 

30 

IF  I IB.LE.OIGO  TO  60 

0008 

CO  50  I«1 ,  IB 

CC09 

50 

IBC-I  IN-ni)4lBC)/I 

0010 

60 

BCCEF-IBC 

0011 

RETURN 

0012 

ENC 

OATE 


68271 


68271 


68271 


-55- 


REFERENCES 

Arnold,  A.  and  Avez,  G.  1967,  Prob limes  ergodiones  de  la  mecanique 
classiaue .  Gauthier-Villars ,  Paris. 

Birkhoff,  G.  1927,  Dynamical  Systems.  American  Mathematical  Society 
Colloquium  Publications,  Volume  IX. 

Bray,  T.  and  Goudas,  C.  L.  1967,  Astron.  c 7.  ,  12 _,  202. 

Brown,  E.  and  Shook,  C.  A.  1933,  Planetary  Theory.  Cambridge  University 
Press,  Cambridge. 

Chai,  A.  and  Kass ,  L.  1966,  AFOSR  Scientific  Report  No.  67-0123. 
Contopoulos,  G.  1963,  Astron.  J.  ,  68^,  11. 

Contopoulos,  G.  1967,  Colloque  sur  les  nouvelles  methodes  de  la 
dynamique  stellaire.  Editions  du  CHRS,  223. 

Deprit,  A.  and  Price,  J.  F.  1965,  Astron.  J.,  70_,  836. 

Deprit,  A.,  Henrard,  J.  and  Rom,  A.  1967,  AAS  Science  and  Technology 

Series,  1J,  Space  Flight  Mechanics  Specialists  Symposium,  85. 

Deprit,  A.  1969,  Celestial  Mechanics ,  1_,  in  print. 

Deprit,  A.  and  Henrard,  J.  1968,  Adv.  in  Astronomy  and  Astrophysics , 

6,  1. 

Deprit,  A.  and  Rabe ,  E.  1969,  Astron.  J.,  7J*_,  in  print. 

Gustavson,  F.  G.  1966,  Astron.  J.,  l\_t  687. 


-56 


Poincare,  H.  1893,  Lee  Rhodes  nouvelles  de  la  mdcanioue  celeste. 
Gauthier-Villars  (Reprinted  by  Dover  Publ.  1957). 

Siegel,  C.  L.  1941,  Ann.  of  Math.,  42,  210. 

Siegel,  C.  L.  1954,  Math.  Ann.,  128,  144. 

Walker,  G.  H.,  and  Ford,  J.  1969,  Manuscript  submitted  for  publication. 
In  private  communication  from  the  Georgia  Institute  of  Technology, 
January  1969. 

Whittaker,  E.  T.  1916,  Pvoo.  Roy.  Soo.  Edinburgh,  37,  97. 


